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Abstract 

In this paper a periodic parameter switching scheme is applied to the Hindmarsh-Rose neuronal system 
to synthesize certain attractors. Results show numerically, via computer graphic simulations, that the 
obtained synthesized attractor belongs to the class of all admissible attractors for the Hindmarsh-Rose 
neuronal system and matches the averaged attractor obtained with the control parameter replaced with 
the averaged switched parameter values. This feature allows us to imagine that living beings are able to 
maintain vital behavior while the control parameter switches so that their dynamical behavior is suitable 
for the given environment. 
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1 Introduction 

Based on the theory of dynamical systems, Hindniarsh and Rose proposed the phenomenological neuron 
model, which may be seen either as a generalization of the Fitzhugh equations or as a simplification of the 
physiologically realistic model proposed by Hodgkin and Huxley [1] [5] . The Hindmarsh-Rose (HR) model 
of neuronal activity is aimed to study the spiking-bursting behavior of the membrane potential observed in 
experiments of a cell in the brain of the pond snail [1 , . The dynamics of a single HR neuron has been studied 
well, and it is illustrated that it can exhibit complex dynamical behavior including periodic and chaotic 
spiking (bursting) motion when certain control parameters of nervous cell models are varied [U [H [3] . 

It is well accepted that the HR neuron model is an alternative candidate for studying the dynamics of 
neuronal systems since it has simpler mathematical forms than Chay's and Hodgkin and Huxley's neuron 
models [HIS]. Therefore, this model has been used to study different aspects of neuronal dynamics such as 
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transitions between different firing regimes [B], relations between block structured dynamics and neuronal 
coding [7] , the effect of noise on neuronal signal transduction |8] and on synchronization J9l , the collective 
dynamics of neuronal networks |10| and the evolution of spiral waves in coupled Hindmarsh-Rose neurons 

m- 

However, an important phenomenon in neuron activity is the transition between different firing patterns. 
From the viewpoint of dynamical systems, it is crucial to investigate the transition mechanism of firing 
patterns for understanding realistic neuronal activities. The local stability and the numerical asymptotic 
analysis of the Hindmarsh-Rose model are used to investigate the firing evolution of a single Hindmarsh-Rose 
neuron |12| . A simple one-dimensional map method has been applied to the HR neurons to convert irregular 
spiking and chaotic bursting to regular beating and periodic bursting [13) . 

Recently, we developed a parameter switching method to synthesize attractors of a class of dynamical 
systems, called hereafter attractors synthesis (AS) algorithm. This method in fact starts from a set of given 
parameter values, and allows us, via periodic or random parameter-switching, to generate any of the set of all 
possible attractors of a class of continuous time dynamical systems of integer order[T ^ .|15 ) . or of fractional 
order [16) . depending linearly on the control parameter. It has been applied successfully to several dynamical 
systems such as Lorenz, Chen, Rossler, Lotka-Volterra, minimal networks, fractional Lii systems and so on. 
Extending this subject, we will synthesize attractors of the HR neuronal model, which can exhibit burst 
dynamics, via the AS method contributing to a better understanding of neuronal firing transition. 

The paper is organized as follows: the HR model will be described in Sec. 2; the AS method will be 
introduced in Sec. 3, and attractors will be synthesized in Sec. 4. Finally, the paper will end in Sec. 5 with 
some comments and conclusions. 

2 Description of Hindmarsh-Rose model 

The dynamics of an isolated HR neuron is governed by the following set of differential equations [T] : 

Xi — bx\ — ax\ + X2 — x^ + I , 

X2= C - dx\ - X2, (1) 
X3^ p[s{xi ~Xi) ~ X3], 

where xi is the membrane potential, X2 is associated with the fast current, Na+, or K+ and 0:3 with the slow 
current, for example, Ca^+. The parameters are [1]: a — I, 6 = 3, c = 1, d = 5, s = 4, xi = — l/2(l + -\/5) ~ 
— 1.6. / and p are the control parameters, and / is a slow parameter while p is a fast parameter. / mimics 
the membrane input current for biological neurons; p controls the speed of variation of the slow variable. 

To find the fixed points we have to solve numerically the following equations: X2 — —hx\ -1-1, x^ — Axi -\- 
6.4, and + 2x1+^X1 + 2^0 which gives the single real solution xi = -0.639, X2 = -1.041, x^ = 3.8442. 
Therefore the system H]) has a single equilibrium point: X* = (-0.639, -1.041, 3.844)^ 

Parameter p is the ratio of time scales between fast dynamics and slow dynamics. Therefore, it controls 
the difference between the slow and the fast dynamics of HR neuron model corresponding to the difference 
between fast fluxes of ions across the membrane and slow ones. Therefore, it is really interesting to investigate 
the dynamics of the HR neuron as the parameter p is changed. However, with / as control parameter, the 
underlying dynamical system does not belong to the class of systems where AS can be applied [see Section 
3]. _ 

The bifurcation diagram of the flrst component xi versus p is shown in Fig. la,b,c for / = 3.4, while 
that for / = 3.5 in Fig. Id. We are interested in the case / = 3.4 because this case represents interesting 
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dynamics, where chaos appears in a narrow range of p with a width of about 3.5 x 10~^. Fig. la gives the 
bifurcation diagram with respect to the control parameter p in the range [0.0001, 0.15]. It is shown that as 
the parameter p is changed, the HR neuron firstly bifurcates from period-doubling to chaos, and then, it 
is stopped via the inverse period-doubling. For a clear vision, the enlargement of Fig. la is shown in Fig. 
lb, which further confirms the above observation. Further enlargement can clearly guide our assessment as 
illustrated in Fig. Ic. 

More details on the dynamics of the HR model can be found in [12] 

3 Attractors synthesis algorithm 

The AS algorithm can be applied to the following general class of continuous-time autonomous and dissipative 
dynamical systems, modeled by the Initial Value Problem (IVP) [15] 

S : fp{x), x{0) = xq, (2) 

where fp is an R"-valued function with a bifurcation parameter p e R, n > 3, and has the expression 

fp{x) ^g(x) +pAx, (3) 

g : M" — > M" is a continuous-time nonlinear function, A is a real constant n x n matrix, xq G M", and 

This class of dynamical systems contains the best-known systems such as Lorenz, Chen, Rossler, Lotka- 
Volterra, minimal networks, fractional Lii systems and so on (see [H]|15j[TB]). 

In [15 it has been shown, via numerical analysis and computer graphic simulations, that the AS algorithm 
allows the synthesis of any attractor of S by parameter switching following some designed rule. 

The AS algorithm can explain what happens with a system when, intentionally or not, the parameter 
value switches quickly through a set of values. Thus, whenp is switched following some designed deterministic 
[15j or random [14! rule, while the system evolves in time, an attractor belonging to the set of attractors is 
generated (synthesized) . 

The AS algorithm can be useful in the cases where some desired control parameter value cannot be 
directly accessed and we want to obtain the corresponding attractor. 
Let us next suppose the following assumptions hold. 

Assumption A.l. Throughout, the existence and uniqueness of solutions of the IVP ©-([S]) are assumed. 

As known, the computer graphic simulations of the numerical integration results of ©-(O can give 
excellent approximations to the orbits within the invariant sets [17j. Thus, the orbits which start near a 
hyperbolic attractor will stay near and they will be shadowed by orbits within the attractor. This happens 
because attractors arise as the limiting behavior of orbits. Therefore, the shadowing property of hyperbolic 
sets [18j enables us to recover long time approximation properties of numerical orbits such as HR's case. 

The AS algorithm consists in using a time varying, or more precisely, a periodically switching parameter, 
according to some design rule. It will be shown, empirically by various experiments, that a desired targeted 
attractor can be duly obtained by the proposed switching scheme. 

Remark 1 The algorithm is robust to some extent: the switching timing and switching parameter values 
both allow flexibility. Therefore, they do not need to be very precisely determined. 

Hereafter the following notations will be used 
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Notation 1 Let us denote by A the set of all global attractors of S, including attractive stable fixed points, 
limit cycles and chaotic (possibly strange) attractors; V cM the ordered set of the corresponding admissible 
values of p and Vn = {pi,P2, ■ ■ ■ ,Pn} C V a finite ordered subset of V containing N different values p, 
which determines the set of attractors An — {^pi : 1 ■ • ■ j -^pn } C »4. 

Assumptions 

A. 2. V is considered, as in most cases, to be composed of a single real interval and all the values of 
'Pn ~ {pi,P2, ■ ■ ■ ,Pn} for which the system behaves stable and/or chaotic is assumed to be accessible. 

n 

A. 3. S is dissipative i.e. ■ f < 0, where xj ■ f = ^ 9/ (xi, X2, . . . , Xn) /dxi (see e.g. [H]). 

Due to the assumed dissipativity, A is non-empt}!!] and it follows naturally that for the considered class 
of systems, a bijection may be defined between the sets V and A. Thus, giving any p ^ V, a. unique global 
attractor is specified, and vice versa. 

Remark 2 Because in this paper computer simulations are used as the major analytical tool, the w -limit 
set (actually, its approximation \20^ ) is considered after neglecting a sufficiently long period of transients. 
Therefore, by attractors (background on the notion of attractor can be found in 1211) it is appropriate to 
understand in this paper the oj-limit set obtained by a numerical method for ODEs with fixed step size h after 
the transients were neglected. 

Let Pn ~ {pijP2t ■ ■ ^Pn}- The AS algorithm relies on the following deterministic time rule applied 
repeatedly on / 

[{mih)pi, {m,2h)p2,...,imNh)pN], (4) 

where the weights are some positive integers. The algorithm acts as follow: in the first time subinterval 
of length mi ft,, p will have the value pi (i.e. the IVP©-® will be integrated mi steps for p = Pi), for the 
next m2 integration steps, p = p2 and so on until the A^-th time subinterval of length vrij^h where p = pn 
and then the algorithm repeats. In order to simplify the notation in (|4]), for a fixed step size h, the scheme 
([4]) will be denoted next 

[mipi, m2P2, . . . ,TOJV PJv]- (5) 

For example, the scheme [pa, 3pi, 2p2\ represents the infinite sequence oip : p^, 3pi, 2p2,P3, 3pi, 2p2, ■ . ■ which 
means that while the considered numerical method integrates ©-([S]), p switches in each m^/i time subin- 
terval. In other words, the numerical method will integrate one step with p = pa, then three times 
with p = pi, then two steps with p — p2 and so on. 

The AS algorithm for the scheme (O is presented in Fig. 2 

Notation 2 Let us denote the synthesized attractor obtained with AS algorithm, via by A* and by Ap* 
the averaged attractor corresponding to 

N 

Pkmk 

P* = ^ ■ (6) 

k=l 

^ Attractor sets can exist only for dissipative systems because shrinking of the volume in phase space for conservative systems 
is ruled out by Liouville theorem 
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Next, the following notion is introduced 

Definition 1 Two attractors are considered to be identical if their underlying a; -limit sets coincide in the 
phase space, the identity being considered from a geometric point of view in the phase space, aided by computer 
graphic analysis. 

This identity in the case of chaotic attractors will be considered only asymptotically since they are fully 
depicted only after an infinite time. 

Proposition 1 The synthesized attractor A* and the averaged attractor Ap* are identical 

N 

Proof. Let us consider some subset Vn- If we denote au — ruk/ is easy to see that p* is a 

1=1 

N N 

convex combination p* — cikPk, since J2 '^k — ^- Therefore p* belongs within the interval (pi, . . . ,pn), 

k=l k=l 

whatever the values pi are chosen. Also, taking into account the bijection between V and A, we are entitled 
to consider that the same convex structure is preserved into A. Therefore, for whatever switched values of p 
in some subset Vn, A* will belong within the set An considered to be ordered by the mentioned bijection 
i.e. Apt S {Ap-^, Apj^). Next, aided by numerical analysis and via computer graphics, it can be showed that 
A* is identical (in the sense defined above) to Ap* .Therefore A* e {Ap-^,Apj^) . m 
Next, we can formulate the main property of the AS algorithm 

Proposition 2 For whatever considered set Vn, the AS algorithm generates an attractor A* which belongs 
to {Ap-^ , Ap^) . 

Remark 3 The AS algorithm is useful especially when we want to obtain some attractor Ap even the un- 
derlying value p cannot be accessible. 

In this case A can be synthesized by choosing a corresponding set {pi, . . . ,pn) ^ p but p ^ Vn = 
{pi, . . . ,Pn}, o-nd a corresponding scheme ([^. 

The initial conditions play an important role since for a specific value p G V there is a single global 
attractor but which could be composed by several local attractors (see e.g. [S] , [H] , [23] j [21] ) • For example, 
for the Lorenz system for the control parameter r — 2.5 there are three local attractors: the origin and two 
symmetrical fixed points Xi,2 (±2, =F 2, 1.5) while for r = 28 there is a single local attractor which is global 
too, the known Lorenz strange attractor. To avoid these possible difficulties, all the computer simulations 
for A* and Ap- for a particular case, start from the same initial conditions. 

In order to see how the AS algorithm works, let us consider that we want to synthesize the attractor Ap. 
Then we must choose N and the set Vn such that p G (pi, . . . ,pn) {p can be equal or not to one of the 
elements pi, ior i = 2, . . . N — 1) . Next, choosing empirically the scheme (O, such that the right-hand side of 
^ gives p, the initial value problem is integrated and finally Ap is obtained. 

To underline the identity between A* and Ap- histograms and Poincare sections besides the phase plots were 
drawn after the transients were neglected. 

4 Synthesis of HR attractors 

Real neurons often display high nonlinearity, which has been shown in many experiments and confirmed by 
numerical simulation of many neuron models, including the HR and Chay neuron models[TJ |4]. Chaos is a 
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universal phenomenon in certain neurons, such as those in the human brain where, often, the information is 
decoded and transduced by means of chaos. It is very important to study the chaos of neuron from different 
aspects. Next, we will focus on the neuronal firing pattern when p is switched. 
Firstly, it is easy to prove the following 

Proposition 3 The HR system (J^l with p considered as control parameter, belongs to the class of systems 
modeled by the IVP 

Proof. Choosing the following substitution 

yi = xi- xi 

and replacing finally yi with Xi, ([IJ) becomes 

x\ — a\x\ + b\x\ + c\x\ + d\X2 + 61X3 + /i + J, (7) 

xi = a2x\ + 62X1 + C2X2 + ^2, 

xz = p{sxi - X3), 



where 



fli = —a, 

hi = b — 3axi, 

ci ~ xi (26 — 3aa;i) , 

di = 1, 

ei = -1, 

/i = -a^l + bxl, 

a2 = -d, 

62 = — 2(ia;i, 

C2 = -1, 

d2 — dxi- 



Thus, the right hand side of ^ can be written following (0) with 

aixl + bixl + cixi + diX2 + £1X3 + fi+ I 
g{x) = ( a2xl + 62X1 + C2X2 + d2 





A- 








^ 






• 








and X = 


X2 
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Because, from a numerical point of view, the mathematical model ([T|) presents a more accessible form 
than ([7]), next we will work with ([T]). 

The HR system verifies the assumptions Al, A2 and also A3 for large parameters range taking into 
account that \/ ■ f — — 3aa:f + 2bxi — 1 — p. Therefore, one can apply the AS algorithm to the HR system 
([T|). For this purpose we considered the Standard Runge-Kutta numerical method for ODEs with the fixed 
step size h — 0.005 and, as stated above, / = 3.4. 

For the sake of the brevity, in what follows, only the most relevant cases are considered. Let us firstly 
consider the case N = 2 and Vn = {0.004,0.01} (see Fig. lb) with nii = m2 = 1. Then, using the scheme 
[lpi,lp2] while integrating the initial value problem ([T]), finally, even the attractors Ao.oo4 and Aq.oi are 
stable limit cycles, the obtained synthesized attractor A* is chaotic (Fig. 3a). As can be seen from Fig. 3a, 
the synthesized attractor A* (blue) and the averaged attractor Ap* (red) for p* = (0.004 + 0.01)/2 = 0.007 
given by ^ coincide. It can be seen that p* € {pi,P2) (Fig. lb) and also A* is situated between the 
corresponding attractors Ap-^ and Ap^ . 

Next, if we use TV = 10 values for p defined as follows: pi — 0.0002 + i x 0.0001 and = 1, for 
i — 1,2,..., 10, the synthesized attractor A* is a stable limit cycle which coincides with the averaged 
attractor Ap. with p* = 0.00075 (Fig. 4). 

An interesting case appears for = 2 and Vn — {0.0082,0.008765} and mi = m2 = 1, where pi, and 
P2 are chosen in a very narrow periodic window (its width is about IE — 4 , see Fig. Ic, Fig. 5b and Fig. 
5c). The obtained attractor. A*, is identical to Ap. for p* — 0.0084825 (Fig. 5 a). We found numerically 
that the fixed point X* is unstable for p — p* since the three eigenvalues are: —6.266,0.179 and 0.02 and, 
as showed in Fig. 6, the Lyapunov spectrum has two negative exponents while the maximum one is positive 
(approximately zero). All these lead us to consider that Ap* (and consequently A*) is a stable limit cycle. 
The apparent difference between the attractors in this case is due to the very small difference between pi 
and p2 ■ 

5 Conclusions 

In summary, we have investigated the synthesis of attractors of the HR neuronal system by means of the AS 
method. It is shown that when we choose the slow parameter p as the control parameter, the HR neuronal 
system belongs to a class of systems, where we can apply the proposed switching method. Consequently, we 
concluded that every attractor can be synthesized by the proposed periodic parameter switching scheme in 
the HR neuron model. It is accepted that neurons can code information by means of firing patterns, which 
depends on the variation of key parameters of neuronal systems. Hence, the results presented in this paper 
may be instructive in understanding the implications of realistic neuronal dynamics. 

Acknowledgements This work was supported by the National Science Foundation of China (Fund 
Nos. 10972001, 10702023 and 10832006). 
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Figure 1: (a) Bifurcation diagram: xi vs p with / = 3.4. (b) Detail of Fig. la. (c) Consecutive detail of Fig. 
lb. (d) Bifurcation diagram: xi vs p with / = 3.5. 
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t = 
repeat 

for k = 1 to N do 
P^Pk 

for i = 1 to mfc do 

integrate Q - © 

t = t + h 

end 

end 
until t > T 

Figure 2: Pseudocode of AS algorithm. 




Figure 3: AS algorithm with scheme [lpi,lp2] with pi — 0.004 and p2 = 0.01; a) The synthesized and 
averaged attractors A* and Ap. respectively, superimposed forp* — 0.007 ; b) Attractor Ao.oo4 ; c) Attractor 
^col- 



li 



A ■, a' with p'=0.00075 



3.S85 - 



3.58- 




Figure 4: AS algorithm for N = 10, p» = 0.0002 + i x 0.0001 and m, = 1, for i = 1,2, . . . ,10; A* and 
Apt, with p* = 0.00075, are superimposed. 
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Adtnaat 








1; a) A* and Ap', with 



Figure 5: AS algorithm for TV 2 and Ptv = {0.0082, 0.008765} and toi = m2 
p* = 0.0084825 , superimposed; b) A0.0082; c) ^0. 008765; d) Superimposed histograms of A* and Ap*; e) 
Superimposed Poincare sections of A* and Ap* for 1.2. 
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Figure 6: Temporal evolution of Lyapunov exponents for p — 0.0084825. 
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